arXivil506.05303v2 [astro-ph.GA] llJul2015 


Mon. Not. R. Astron. Soc. 000, 1-14 (2015) Printed 14 July 2015 (MN lATgX style file v2.2) 

Dynamical Constraints on the Origin of Multiple Stellar 
Populations in Globular Clusters 

P. Khalaj* and H. Baumgardt 

School of Mathematics and Physics, University of Queensland, St. Lueia, QLD 4072, Australia 


Accepted xxxx. Received xxxx; in original form xxxx 


ABSTRACT 

We have carried out a large grid of Wbody simulations in order to investigate if mass- 
loss as a result of primordial gas expulsion can be responsible for the large fraction of 
second generation stars in globular clusters (GCs) with multiple stellar populations 
(MSPs). Our clusters start with two stellar populations in which 10% of all stars are 
second generation stars. We simulate clusters with different initial masses, different 
ratios of the half-mass radius of first to second generation stars, different primordial 
gas fractions and Galactic tidal fields with varying strength. We then let our clusters 
undergo primordial gas-loss and obtain their final properties such as mass, half-mass 
radius and the fraction of second generation stars. Using our Wbody grid we then 
perform a Monte Carlo analysis to constrain the initial masses, radii and required 
gas expulsion time-scales of GCs with MSPs. Our results can explain the present-day 
properties of GCs only if (1) a substantial amount of gas was present in the clusters 
after the formation of second generation stars and (2) gas expulsion time-scales were 
extremely short (< 10^ yr). Such short gas expulsion time-scales are in agreement with 
recent predictions that dark remnants have ejected the primordial gas from globular 
clusters, and pose a potential problem for the AGB scenario. In addition, our results 
predict a strong anti-correlation between the number ratio of second-generation stars 
in GCs and the present-day mass of GCs. So far, the observational data show only a 
significantly weaker anti-correlation, if any at all. 

Key words: globular clusters: general - stars: chemically peculiar - stars: formation 
- stars: kinematics and dynamics - methods: numerical 


1 INTRODUCTION 

It is generally assumed that all stars in star clusters are born 
in close proximity to each other, and in well-mixed molec¬ 
ular clouds by a rapid star formation process and therefore 
have similar ages and metallicities (Lada & Lada 2003). As 
a result, star clusters should only host a single population 
of stars, i.e. they are bona fide single stellar population 
systems. However, recent observations of GCs show a statis¬ 
tically significant star-to-star variation in the abundance of 
light elements, such as Na, O, Mg or A1 (e.g., Carretta et al. 
2009; Gratton et al. 2013). These abundance anomalies of 
light elements are not associated with any spread in the iron 
abundance for the majority of GCs except for a few cases 
such as LJ Cen (Gratton, Sneden & Carretta 2004). Such 
massive GCs are thought to be different than normal GCs 
and have a different origin, for example being the remnant 
of a disrupted dwarf galaxy (Meza et al. 2005). 

In addition to the abundance anomalies mentioned 
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above, the colour-magnitude diagrams of some GCs split 
into two or more evolutionary sequences (e.g., u Cen Rey et 
al. 2004 and Bedin et al. 2004; NGC 2808 D’Antona & Caloi 
2004, Piotto et al. 2007 and Milone et al. 2012; NGC 1851 
Milone et al. 2008; 47 Tuc Milone et al. 2012; NGC 6397 
Milone et al. 2012; M22 Marino et al. 2012; GCs in Fornax 
D’Antona et al. 2013). 

These findings are indicative of self-enrichment in GCs 
and suggest that star clusters are comprised of at least two 
stellar populations, in direct contradiction to the conven¬ 
tional star formation scenario described earlier. We refer 
to these two populations as first generation (FG) and sec¬ 
ond generation (SG) stars for convenience and to be con¬ 
sistent with previous studies such as Decressin, Baumgardt 
& Kroupa (2008). In our terminology FG and SG stars cor¬ 
respond to stars with normal (or primordial) and enriched 
chemical compositions respectively. 

The observations show that the number ratio of SG to 
FG stars, N 2 /Ni, is around unity although with considerable 
spread (D’Antona & Caloi 2008). Further studies on MSPs 
have shown that they exhibit different spatial (u Cen Bellini 
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et al. 2009; several GCs Lardo et al. 2011; 47 Tuc Natal et 
al. 2011; M15 Larsen et al. 2015) and dynamical signatures 
(47 Tuc Richer et al. 2013 and Kucinskas, Dobrovolskas & 
Bonifacio 2014). However, Dalessandro et al. (2014) found 
that the different populations in NGC 6362 share the same 
radial distribution which is the first evidence of fully spa¬ 
tially mixed MSPs ever observed in a GG. 

Several scenarios have been proposed to address the 
origin of multiple stellar generations in GGs (Decressin et 
al. 2007; Decressin, Gharbonnel & Meynet 2007; de Mink 
et al. 2009; Renzini 2008; D’Ercole et al. 2010; Gonroy & 
Sprgel 2011; Ventura et al. 2001; Valcarce V Gatelan 2011 
and Bastian et al. 2013) among which the four main scenar¬ 
ios are: (1) fast rotating massive stars (FRMS) (20 —120 M©; 
Prantzos V Gharbonnel 2006; Maeder V Meynet 2006; De¬ 
cressin et al. 2007; Decressin, Gharbonnel V Meynet 2007); 
(2) asymptotic giant branch (AGB) stars (4 —OM©; Ventura 
et al. 2001; D’Ercole et al. 2008, 2010; Ventura V D’Antona 
2011); (3) massive (10 — 100M©) and intermediate mass 
(4— 10 M©) binaries (de Mink et al. 2009) and (4) early disc 
accretion in low-mass pre-main-sequence stars (enriched gas 
comes from stars with M > 10M0)(Bastian et al. 2013). 

In the ERMS scenario, stars that spin at a rate close to 
their critical break-up speed can lose extensive amounts of 
mass via stellar winds which are slow enough to be retained 
in the gravitational potential well of the cluster and form cir- 
cumstellar discs out of which SG stars will be born. In the 
AGB and massive binary scenarios, these slow winds come 
from the envelopes of evolving intermediate-mass stars in 
their AGB phase and numerous massive interacting binaries 
in the core of clusters respectively. SG stars in the FRMS 
scenario need to form in a very short time-scale, t < 8.8 Myr 
Krause et al. (2013), before the burst of the first supernova 
(SN), since SN winds can destroy circumstellar discs formed 
around fast-rotating EG stars and interrupt the star forma¬ 
tion (Decressin et al. 2007). In the AGB scenario, on the 
other hand, the formation of SG stars is triggered after all 
SNe have gone off and the cluster has been cleared of SN II 
ejecta (t > 28 Myr) (D’Ercole et al. 2008), otherwise AGB 
ejecta will be polluted by SN ejecta which will cause the iron 
abundance of SG stars to differ from that of EG stars, in 
contradiction with observations. Bastian et al. (2013) pro¬ 
posed a model in which GGs do not need to go through 
different instances of star formation to produce chemically 
peculiar stars. According to this model, interacting massive 
binaries (M > 10 M©) supply the intra-cluster medium with 
enriched material which will be accreted by low-mass stars 
(M < 2M©) while they are still in their pre-main-sequence 
phase. The main difference between this model and other 
models is that stars with different chemical abundances be¬ 
long to the same generation of stars. The main caveat of 
this model is that the circumstellar discs around accreting 
low-mass stars need to survive for 5 to 10 Myr which is a 
questionable assumption in GGs with a denser core (Bas¬ 
tian et al. 2013). 

In a recent study, Bastian, Gabrera-Ziri V Salaris (2015) 
tested the yields of all the proposed models in the litera¬ 
ture for their consistency with observations and concluded 
that none of the models is able to explain the observed He 
abundance of clusters. As a result the origin of abundance 
anomalies in GGs is still a matter of debate. 

In addition to this discrepancy between the theoreti¬ 


cal yields and observations, the ejecta in the FRMS and 
AGB scenarios are not enough to form a large population of 
SG stars and explain the roughly equal number of EG and 
SG stars found in observations (D’Antona V Galoi 2008). 
Assuming a canonical Kroupa (2001) initial mass function 
(IMF), the total mass which is lost by all FRMS and AGB 
stars, constitutes between ~ 4% to 9% of the initial mass 
of all EG stars (de Mink et al. 2009). If we assume that the 
gas which is lost by this mechanism entirely turns into SG 
stars, i.e. star formation efficiency is 100%, the number ratio 
of SG stars to EG stars is expected to be at almost ~ 10%. 
This issue, which is referred to as the mass-budget problem, 
does not exist for the massive and intermediate-mass bina¬ 
ries as they provide more ejecta than AGB and fast-rotating 
massive stars combined. In addition their ejecta are further 
mixed with an approximately equal amount of pristine gas 
which doubles the mass of the available gas for star forma¬ 
tion. As a result there is a substantial amount of polluted 
gas to form a large number of SG stars (de Mink et al. 2009). 

To address the mass-budget problem two solutions have 
been proposed: either a cluster must have been at least 10-20 
times more massive and have undergone significant mass- 
loss (~ 90%) or the cluster IMF must have been strongly 
top-heavy, i.e. it initially had many more massive stars than 
predicted by a canonical IMF. 

There is observational evidence against both of these 
solutions. First, observations of GGs in a number of dwarf 
galaxies show a high ratio of metal-poor GGs to field stars 
which cannot be explained if star clusters were initially 
10 times more massive and underwent significant mass-loss 
(Larsen, Strader & Bordie 2012; Larsen et al. 2014). Second, 
Dabringhausen, Kroupa & Baumgardt (2009) found that a 
top-heavy IMF will lead to high mass-to-light (M/L) ra¬ 
tios in old stellar systems (t = 12 Gyr) such as UGDs and 
GGs. This is a serious issue for the AGB scenario, as the 
polluters in the AGB scenario evolve into white dwarfs and 
the retention factor of white dwarfs is very high compared 
to the FRMS scenario in which polluters evolve into black 
holes or neutron stars, many of which will leave the cluster. 
In the AGB scenario, depending on whether SG stars form 
as a distinct generation or they are only contaminated by 
the processed gas during formation, one needs a high-mass 
slope of a = —1.15 and a = —1.95 respectively^ to pro¬ 
vide enough ejecta from AGB stars to form low-mass stars 
(Scenarios I and H of Prantzos & Gharbonnel 2006). Ac¬ 
cording to Fig. 2 of Dabringhausen, Kroupa & Baumgardt 
(2009), this will translate into a normalised M/L ratio of 
5.3 and 4.2M©L0^ if the retention factor of SN remnants 
is 0 and 6.5 and 4.3 M©Lq^ if it is 20%. The average ob¬ 
served M/L ratio of GGs in the Milky Way, its satellites and 
M31 is less than 2.0 M© Lq^ (McLaughlin & van der Marel 
2005; Strader, Galdwell & Seth 2011). If one considers the 
biases that exist in the derivation of masses from integrated 
light of GGs, the observed M/L ratios can be explained by 
a canonical IMF (Shanahan Sz Gieles 2015). 

The focus of the present paper is to study the effect of 
significant mass-loss on the dynamical evolution of star clus¬ 
ters with MSPs. The dynamical effects of a top-heavy IMF 
on GGs with MSPs can be further examined in a future pa- 


^ The high-mass slope of a Kroupa (2001) IMF is o: = —2.35 
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per. Using A/’-body simulations we determine the required 
initial conditions under which the final number of SG to 
FG stars match the observations. We then perform a Monte 
Garlo (MG) analysis to compare the outcome of our simula¬ 
tions with observations and determine whether the signifi¬ 
cant mass-loss scenario is able to explain the observed mass 
and half-mass radius distributions of GGs and ultimately be 
the reason for the observed abundance anomalies. 

The present paper is structured as follows. In the next 
section we briefly review the mass-loss mechanisms which 
can affect the dynamical evolution of GGs. We discuss the 
details of our Wbody simulations and the procedures we 
have followed to create our grid of runs in Section 3. In 
Section 4 we present our results and then compare the out¬ 
come of simulations with observations using a MG analysis 
in Section 5. We finally conclude our work in Section 6. 


2 MASS-LOSS MECHANISMS 

There are several mechanisms through which a star clus¬ 
ter can lose mass: stellar evolution induced mass-loss, two- 
body relaxation, external tidal shocks and primordial gas 
loss. These mechanisms are discussed further below. 


2.1 Stellar evolution induced mass-loss 

To lose a large amount of mass via stellar evolution, clusters 
need to be very extended and in a strong tidal field where 
the ratio of the tidal radius to the half-mass radius, 
is small. The effect of stellar evolution induced mass-loss in 
the AGB scenario has been studied by D’Ercole et al. (2008) 
using a series of A-body simulations. In their models, they 
use King (1966) models with Mfg = 10^ M©, vt = 200 pc 
and Wo = 7.0, c = log(rt/rc) = 1.50, where SG stars are 
highly concentrated in the innermost regions of the clusters 
with a half-mass radius one-tenth of that of the initial FG 
stars. Such initial parameters correspond to very extended 
clusters (rtjrh — 8.75 and rn = 23 pc) in which SG and 
FG stars are dynamically decoupled from each other. As 
a result FG stars can readily expand their orbits and be 
stripped by the Galactic tidal field in response to a small 
amount of mass-loss, leaving a sub-cluster of SG stars in the 
center of the initial cluster. D’Ercole et al. (2008) also study 
other models with different initial truncation radii and con¬ 
centrations which underfill their Roche lobes, but the only 
models that match the observed number ratio of SG to EG 
stars, i.e. N 2 /N 1 ^ 1.0, are the very extended and tidally 
filling clusters. This is not in agreement with the observa¬ 
tion of young massive star clusters and today’s properties 
of GGs as they have typical half-mass radii of around 1.0 pc 
(Portegies Zwart, McMillan & Gieles 2010) and 5.0 pc re¬ 
spectively (Harris 1996). Unless the condition under which 
GGs have formed were significantly different, the stellar evo¬ 
lution induced mass-loss cannot lead to significant mass-loss. 
In addition, it is unclear if a sample of clusters with these 
initial conditions can explain the observed distribution of 
SG number ratios in GGs. 


2.2 Two-body relaxation 

Baumgardt & Makino (2003) studied the effect of two-body 
relaxation as well as stellar evolution on the dynamical evo¬ 
lution of star clusters in external tidal fields through A-body 
simulations. They assumed different Galactic orbits, stellar 
density profiles and particle numbers for star clusters and 
derive the following formula for the lifetime of a star cluster 


Tdiss 

Myr 


A 

X 

Kg 

( ] 

ln(0.02iV)_ 


i^220kms“V 


where A is the number of particles, Vb and Rg are the 
Galactic circular velocity and distance of the cluster and 
e is the eccentricity of the cluster orbit, x and [3 are two 
parameters whose values depend on the initial concentration 
of the cluster and for King Wo = 7.0 they are equal to 0.82 
and 1.03 respectively. 

Eor A = 10®, e = 0.5, Vq = 220kms“^ and Rq = 
8.5 kpc, Tdiss will be about 55Gyr which shows that two- 
body relaxation is a slow process for massive GGs and is 
not efficient in reducing the mass of GGs by 90% over one 
Hubble time. As a result this process cannot be the origin 
of significant mass-loss in star clusters and we will omit this 
process in our A-body simulations. 

Eor a cluster whose initial number ratio of SG to EG 
stars is ~ 10% and SG stars are more concentrated than 
EG stars, two-body relaxation causes different stellar popu¬ 
lations to fully mix in about 2 elapsed half-mass relaxation 
times (Decressin, Baumgardt & Kroupa 2008). Using Eq. 1 
of their paper, a cluster with M = 10® M© and = 3 pc, 
has a mixing time of approximately 2 Gyr. This implies that 
any significant mass-loss scenario proposed to explain the 
origin of MSPs must have a shorter time-scale, since after 
the mixing has occurred the number ratio of SG to EG stars 
will not change due to further mass loss. 


2.3 External tidal shocks 

External tidal shocks such as encounters with giant molec¬ 
ular clouds (GMGs) are able to disrupt open clusters M < 
10^ M0 via a single encounter on time-scales of about ~ 
2.0 Gyr (Widen 1985; Gieles et al. 2006). Gieles et al. (2006) 
derived the following formula for the disruption time of star 
clusters 


Tdis = 2.0 


5.1 M| pc~® 
EnPn 


Me 


104 M© 


0.61 

Gyr 


where Sn and pn are the individual surface and global den¬ 
sity of the GMGs, equal to 17OM0pc“^ and O.O3M0pc“® 
in the solar neighbourhood (Solomon et al. 1987). Eor a GG 
with Me = 10® M©, this formula gives a disruption time of 
almost ~ 33 Gyr. In denser environments such as the center 
of M51, Pn is 10 times higher (Gieles et al. 2006) which short¬ 
ens the disruption time by an order of magnitude, but this 
is still larger than the mixing time of MSPs (~ 2Gyr), as 
discussed in Section 2.2. In addition encounters with GMGs 
are stochastic by nature. Hence they cannot be responsi¬ 
ble for significant mass-loss in all clusters and their effect is 
insignificant over short time-scales. 
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2.4 Primordial gas loss 

If the star formation efficiency is less than 100%, this pro¬ 
cess will happen to every cluster of any size or mass since 
it has an intrinsic origin. Any gas loss in GCs will be ac¬ 
companied by loss of stars, especially when the gas loss is 
impulsive (Baumgardt & Kroupa 2007). There are a number 
of different sources which can inject enough energy into the 
intra-cluster medium to entirely unbind the primordial gas. 
Examples are stellar winds, SN explosions (Decressin et al. 
2010) and black holes (Krause et al. 2012, 2013). 

With all other mass-loss scenarios excluded or shown 
to be ineffective on short time-scales, primordial gas loss 
remains as the only plausible and universal mechanism via 
which GCs can lose a significant amount of mass over a few 
Myrs and in our A’-body simulations we only deal with such 
primordial gas loss as discussed in the next section. As men¬ 
tioned above accretion onto dark remnants is one candidate 
for a mechanism which can cause such a primordial gas loss. 
The setup of our model clusters and our analysis, though 
consistent with the dark remnant scenario (Section 5), is 
not limited to this scenario and in principle can be applied 
to any other physical process that has a similar effect on 
GCs. 


3 iV-BODY SIMULATIONS 


We set up clusters consisting of 3 components: EG stars 
(~ 90% of total stellar mass), SG stars (~ 10%) and a gas 
cloud whose mass is a free parameter in the simulations. 
The initial number ratio of SG to EG stars N 2 /N 1 is fixed 
at 0.1. We do not directly simulate the gas particles but 
only calculate the force that the gas cloud exerts on each 
star. The initial density profiles of the different components 
are given by Plummer (1911) models with different masses 
and Plummer radii. We have used Plummer models since 
they are easy to work with and it is also possible to validate 
the outcome of the simulations using analytical methods. We 
do not expect that other initial density distributions such as 
King (1966) models will affect the final results significantly, 
as the exact details of any initial density distribution will 
be quickly wiped out by violent relaxation as a result of 
significant mass-loss in the simulated star clusters. 

The central gas cloud and SG stars have the same degree 
of concentration with respect to EG stars in our simulations, 
i.e. a 2 /ai = %/ai G {0.1, 0.2} where ai, a 2 and ag are the 
Plummer scale radii of EG stars, SG stars and the gas cloud 
respectively. Table 1 summarizes the initial conditions of our 
simulations. 

The time-scale of our simulations is short compared to 
the clusters relaxation times so the effect of stellar evolu¬ 
tion, mass segregation, etc. can be neglected. As a result 
all particles in our simulations have equal masses which are 
constant throughout the whole simulation. All time-scales 
in our simulations are expressed in terms of the the initial 
crossing time of the cluster which is defined to be: 


Ter = 


2r/r 
(7 V 


( 1 ) 


Where ru is the initial half-mass radius of all stars (which 
is approximately equal to 1.20 times the Plummer radius of 
EG stars for the values adopted in Table 1, i.e. ru ~ 1.20ai) 


and av is the initial velocity dispersion of stars calculated 
from the virial theorem in the presence of gas. 

We start with clusters which are initially in virial equi¬ 
librium and then remove the gas according to the following 
equation 

( 2 ) 

[MgiO) t<t0 


where t is the simulation time, to is the amount of time that 
we wait before removing the gas^ and is set to be equal to 
5 crossing times in all simulations and r = Texp/Tcr is the 
ratio of gas expulsion time-scale to the initial crossing time. 
We change r in the range 10“^ to 10^ on a logarithmic scale, 
corresponding to instantaneous and adiabatic gas expulsion 
respectively. 

The initial mass of gas Mgas(O) is parameterized by a 
parameter 77 which is the ratio of the initial mass of gas 
divided by the initial mass of EG stars, i.e. 


MAP) 

' Mi(0) 


( 3 ) 


where 77 varies from 0 to 2.00 in steps of 0.02 in our simula¬ 
tions. 

All simulated clusters are in a Galactic tidal field which 
is modelled using the near-held approximation (Aarseth, Lin 
& Palmer 1993) and implemented by writing the equations 
of motions of stars in a right-handed rotating coordinate 
system whose origin is initially centered on the cluster and 
X and y axes point toward the Galactic anti-center and the 
direction of orbital motion of the star cluster respectively, 
assuming that the star cluster is moving in the x — y plane. 
The equation of motion for a star in such a coordinate sys¬ 
tem is given by 

(total) = K(stars) + ri(gas) — 2ft x K + ^“^{Sxiex — ZiQz) 

( 4 ) 

Where (stars) + (gas) is the acceleration of each star due 
to the total gravitational force of other stars and the gas 
cloud which are calculated using the following equations 


j = N 

K (stars) = 


Grrij 


Ir. - r^P + 


3/2 


{Vj - n) (5) 


ri(gas) 


GMgjt) 


rf + 


3/2 * 


( 6 ) 


where e is the softening parameter that we have introduced 
in our simulations and it is equal to the minimum distance 
between stars in the central region of the cluster. The third 
and forth terms on the right-hand side of Eq. (4) are the 
Goriolis and centrifugal force combined with the tidal forces 
respectively and Q. = is the angular velocity of the 

cluster around the Galactic center. 

In our simulations the strength of tidal field is param¬ 
eterized by the ratio of the tidal radius of each cluster to 


^ The reason that we don’t remove gas at t = 0 is that we want 
to measure the dynamical properties of the simulated clusters in 
the first few crossing times when the cluster is still in equilibrium. 
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Table 1. Initial conditions of the simulations. 


Parameter 

Symbol 

Value 

Fixed Parameters 



Total number of stars 

N 

20480 

Initial number of SG to FG stars 

N 2 

0.1 

Ratio of the Plummer scale radius of the gas cloud to SG stars 

ttg 

1.0 

a2 


Variable Parameters 



Ratio of the Plummer scale radius of SG stars to FG stars 

A=^ 

{0.1,0.2} 


ai 

Ratio of the Initial mass of the gas cloud to the mass of FG stars 

II 

0.0 < 7? < 2.0 

Ratio of the gas expulsion time-scale to the initial crossing time 

^ Texp 

10-2 < T < 10"^ 

Ratio of the initial tidal radius to the Plummer radius of FG stars 

n 

{5,10,15, 20, 25, 30, 35,40, 00} 


ai 


the Plummer radius of FG stars rtjai. We vary rt/ai from 
5 to 40 in steps of 5, where a large value of rt/ai means 
a weak tidal field. We also did one set of simulations for 
rt/ai = oo (Q = 0 ) which corresponds to isolated clusters. 
rt is related to the total cluster mass M*(t) + Mgas(t) and 
Q via the following equation (Giersz & Heggie 1997): 

= (V) 

As a result, all clusters in our grid can be modelled by 
only 4 parameters 77 , n/ai, r and A = a 2 /ai (see Table 1). 
Our simulations can be thought of as a generalized version 
of the Baumgardt & Kroupa (2007) models who considered 
the effect of gas expulsion on a single stellar population. 

All simulations in our grid are run for 555 initial cross¬ 
ing times which was found to be enough for clusters with 
T < 10 ^ to end up in a quasi-equilibrium state after which 
we determined the mass, half-mass radius and number ra¬ 
tio of SG stars. For r > 10^ the gas expulsion is adiabatic 
which only affects the cluster in long term and is dealt with 
in Section 5. 

We ran all the simulations on the Green II GPU super¬ 
computer at Swinburne University of Technology. 


4 RESULTS 

We record the properties of our model clusters throughout 
the simulation. In particular, we find unbound stars using 
an iterative algorithm as described below: 

(i) Find the coordinate of the cluster density center using 
the von Hoerner (1963) method and the unbiased density 
estimator of Gasertano & Hut (1985) for the 10 th nearest 
neighbours of each star, i.e. j = 10 in equation 11.2 of Gaser¬ 
tano & Hut (1985). 


(ii) Using Eq. (9), calculate the instantaneous tidal radius 
of the cluster rt{t) as a function of the remaining cluster 
mass. 

(iii) Find the stars whose distances are larger than the 
tidal radius calculated in the previous step and mark them 
as unbound stars. For isolated clusters (n/ai = 00 ), use the 
total energy of each star as the selection criterion. 

(iv) Subtract the mass of unbound stars from the total 
mass of the cluster. 

(v) Repeat the previous steps until all bound stars reside 
within the tidal radius or all stars are designated as unbound 
stars (i.e. total disruption). 

Using the above algorithm we calculate the mass-loss, 
number ratio of SG stars and expansion factor for the model 
clusters at each instant of the simulation. The outcome of 
our simulations is shown in Fig. 1 to 3 which depict the mass- 
loss AM/Mi, number ratio of SG stars ^ 2/(^1 + N 2 ) and 
logarithm of expansion factors log^Q {rhf /rhi) as a function 
of gas fraction 77 and the ratio of the gas expulsion time-scale 
to the initial crossing time r for different tidal radii ratios 
{rt/ai = 5,10,15, 20, 25,30,35,40, 00 ) and A = (0.1,0.2). In 
these figures, the region enclosed in solid lines corresponds 
to 90 ±5% mass-loss and a value of 50± 10% for the fraction 
of SG stars. The white filled area in the lower right corner 
of each plot is the total disruption zone in which all clusters 
will be totally destroyed as a result of significant mass-loss. 
The region between the black solid lines represents the set of 
initial conditions which match the observations. One can see 
that the width of this region increases in stronger tidal fields 
and for higher concentrations of SG stars (e.g. A = 0.1). 

The trend that we see in these figures can be explained 
as follows: First, the loss of gravitational potential is greater 
for clusters with a higher gas fraction. Second, very short 
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Figure 1. Mass-loss AM/Mi as a function of gas fraction 77 and gas expulsion time-scale r for different strengths of the tidal field 
(pt/cLi = 5,10,15, 20, 25, 30, 35, 40, 00 ) and concentration of SG stars A. The top and the bottom panels correspond to A = 0.1 and 
A = 0.2 respectively. The region enclosed by solid black lines corresponds to 90 ± 5% mass-loss. The white filled area in the lower right 
corner of each plot shows the region where all clusters are totally destroyed. 
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Figure 2. Same as Fig. 1 but for the number ratio of SG stars N 2 /{Ni + N 2 ). The region enclosed by solid black lines corresponds to 
a number ratio of 50 ± 10% for SG stars. 
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rj 



Figure 3. Same as Fig. 1 but for the logarithm of the expansion factor of the cluster log^^Q {'^hf: for A = 0.1 (top) and A = 0.2 
(bottom). 
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Table 2. Range of the initial parameters used in the MC simu¬ 
lations. 

Parameter 

Range 

Steps 


[5.4, 6.5] 

0.05 


[0.7, 1.5] 

0.05 

log(^) 

[3.2, 5.2] 

0.05 

log(^)“ 

[-0.25, 0.5] 

0.05 

‘^log(Texp/yr) = <^log(M*/M0) = 

[0.25, 0.75] 

0.25 

^^og(rh/ pc) 

[0.15, 0.45] 

0.15 

^ This parameter is only relevant 

in models 

with mass- 


independent radii. 


gas expulsion times-scales do not allow loosely-bound stars 
(mainly FG stars) to compensate for the loss of gravita¬ 
tional potential energy and go into an equilibrium so they 
leave the cluster after a few crossing times, whereas in mod¬ 
els with longer gas expulsion time-scales stars have enough 
time to gradually expand their orbits and remain in a quasi¬ 
equilibrium state without crossing the tidal radius and es¬ 
caping from the cluster. Third, FG stars have a lower con¬ 
centration than SG stars and when the gas expulsion is in¬ 
stantaneous, the cluster preferentially loses more FG stars, 
while in the adiabatic case, many of FG stars will be re¬ 
tained in the cluster. Fourth, clusters that lose a substantial 
number of stars have smaller tidal radii and have shrunk 
in size, hence their expansion factor is less than 1.0 and de¬ 
creases with mass-loss. As a result, the mass-loss is expected 
to be much more extreme for clusters with higher gas frac¬ 
tions and shorter time-scales and such clusters must show a 
higher number ratio of SG stars and relatively lower expan¬ 
sion factors. In addition, there should be a region in the pa¬ 
rameter space in which all clusters will be totally disrupted. 
The outcome of our A’-body simulations are consistent with 
Decressin et al. (2010) who analysed the A-body models of 
Baumgardt & Kroupa (2007). 

As it is inferred from Fig. 1 to 3, the initial conditions 
which meet the observational criteria occupy a very narrow 
strip in the parameter space. In addition, this region is very 
close to the total disruption zone, meaning that if we slightly 
change the initial conditions, we will either end up in the 
total disruption zone (SG fraction ~ 100%) or the region 
which is far from the observed clusters (SG fraction ~ 10%). 
We did a MG analysis as explained in Section 5 in order to 
find the physical initial conditions of GGs in terms of cluster 
mass, half-mass radius and gas expulsion time-scale. 


5 MONTE CARLO ANALYSIS AND 

COMPARISON WITH OBSERVATIONS 

In this section we describe the details of our MG analysis on 
the initial conditions of star clusters. We make different sets 


of initial conditions for star clusters and we feed these initial 
conditions into our grid to hnd the hnal conditions and com¬ 
pare them with observations of Galactic GGs. We change the 
initial distribution until the best match with observations is 
found. We have performed our MG analysis for A = 0.1 
and 0.2 separately. We adopt a log-normal distribution for 
the initial distribution of the cluster stellar mass (Parmen- 
tier & Gilmore 2007, 2008) parameterized by a mean value 
and standard deviation of log(M*/MQ) and cr\og{M^/MQ) 
respectively. We assume similar normal and log-normal dis¬ 
tributions for the gas fraction and gas expulsion time-scale 
with mean values of fj and log (Tgxp/ yr) and standard devi¬ 
ations equal to cr(^) and /yr)* assume that the 

initial half-mass radii and the initial masses of GGs are re¬ 
lated via the following initial mass-radius relation derived 
by Gieles et al. (2010) 

log (^) = -3.5650+ 0.615log (^) (8) 

For comparison, we have also done one set of MG simu¬ 
lations by relaxing the mass-dependent constraint on radii 
and replacing it with a log-normal distribution to see how 
it affects the final results. 

Tidal radii of GGs depend on the environment in which 
they form which is unknown. Possible choices are (1) GGs 
have formed in an environment similar to the present-day 
Milky Way, when most of its mass was already in place or 
(2) they formed in satellite galaxies of the Milky Way with 
many of them being disrupted and merged with the Milky 
Way (Prieto & Gnedin 2008) and some survived like the 
LMG and Fornax dwarf galaxy. For the first case we assume 
that our clusters are in a Galactic field with a constant cir¬ 
cular velocity of Vb = 220kms“^ at a distance Rg from the 
Galactic center. The tidal radius can then be determined us¬ 
ing the following equation (Eq. 1 of Baumgardt & Makino 
2003) 


n 



1/3 


r: 


2/3 


(9) 


In our MG analysis we have considered three cases of 
Rg = 2, 4 and 8.5 kpc (solar neighbourhood) corresponding 
to strong, moderate and weak tidal helds in the present-day 
Milky Way. We will refer to these cases as MD-RG2.0, MD- 
RG4.0 and MD-RG8.5 for mass-dependent radii, as given by 
Eq. (8), and MI-RG2.0, MI-RG4.0 and MI-RG8.5 for mass- 
independent radii. 

If GGs formed in dwarf galaxies, they would have tidal 
radii which are comparable to the Rg = 8.5 kpc case. One 
can take LMG with Rg = 4 kpc, Vb = 70kms“^ (Alves & 
Nelson 2000) or Eornax with Rg = 0.5 kpc, Vb = 10kms“^ 
(Strigari et al. 2006) as an example, where Rg refers to 
the radius of the galaxy and Vb is the circular velocity 
at Rg- Using Eq. (9), the tidal radius of a cluster with 
a mass of M = 10® M© in such galaxies is about 190 
and 175 pc respectively, close to the value of 150 pc for 
Rg = 8.5 kpc, Vb = 220 kms“^. As a result our three choices 
of the tidal strength are sufficient to represent different en¬ 
vironments in which GGs might have formed. 

Given the tidal radius and the initial half-mass radius 
of each cluster, the ratio of rt/vh and consequently rt/ai can 
be calculated. As a result all the required initial conditions 
(?7,T, rt/ai) to identify our model clusters in the A-body 
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Table 3. Outcome of the MC simulations. We have simulated 1000 clusters for each set of the initial parameters. Each row shows the 
best-fitting model in terms of the D parameter which is a measure of the goodness of fit and defined by Eq. (10). All models in this table 
have sample mean values of 50 ± 5% for the fraction of SG stars, r refers to the value of the anti-correlation between the cluster mass 
and the fraction of SG stars. The values reported in this table are the means of 20 samples with different random seed numbers. The 
best-fitting models show statistical fluctuations within ±0.1 for log (M*/ Mq), fj and log (Texp/ yr) which can be inferred as an error-bar 
on the best-fitting parameters. 


Tidal Field A log (^^^^±(7 r] ± a log ± a log (^^^±(7 r D {xlO 


mass-dependent radii 

MD-RG2.0 

0.1 

5.95 ± 0.75 

1.25 ±0.50 

4.45 ± 0.50 

- 

-0.50 

1.34 

MD-RG2.0 

0.2 

6.05 ± 0.50 

1.30 ±0.25 

3.95 ± 0.75 

- 

-0.72 

1.86 

MD-RG4.0 

0.1 

5.60 ±0.50 

1.30 ±0.50 

4.30 ±0.50 

- 

-0.70 

1.33 

MD-RG4.0 

0.2 

5.65 ± 0.50 

1.15 ±0.25 

3.60 ±0.75 

- 

-0.72 

1.61 

MD-RG8.5 

0.1 

5.65 ± 0.50 

1.10 ±0.50 

4.20 ±0.25 

- 

-0.72 

1.34 

MD-RG8.5 

0.2 

5.60 ±0.50 

1.20 ±0.25 

4.15 ±0.50 

- 

-0.70 

1.58 

mass-independent radii 

MI-RG2.0 

0.1 

6.10 ±0.25 

0.95 ± 0.25 

3.90 ±0.25 

0.10 ±0.30 

-0.85 

0.87 

MI-RG2.0 

0.2 

6.15 ±0.25 

1.00 ±0.25 

3.40 ± 0.25 

0.05 ±0.30 

-0.86 

0.99 

MI-RG4.0 

0.1 

6.15 ±0.25 

0.95 ± 0.25 

3.55 ± 0.25 

0.10 ±0.30 

-0.83 

0.90 

MI-RG4.0 

0.2 

6.20 ±0.25 

1.05 ±0.25 

3.40 ± 0.25 

-0.05 ± 0.30 

-0.87 

0.99 

MI-RG8.5 

0.1 

6.05 ± 0.25 

0.95 ± 0.25 

3.65 ± 0.25 

0.05 ±0.30 

-0.84 

0.87 

MI-RG8.5 

0.2 

6.10 ±0.25 

1.05 ±0.25 

3.40 ± 0.25 

-0.05 ± 0.30 

-0.86 

0.97 


grid will be uniquely determined and we can find the final 
properties of the clusters by interpolation between the val¬ 
ues of the grid. In order to interpolate in our grid we need to 
assume that mass-loss, fraction of SG stars and expansion 
factor in the total disruption zone (white area in Fig. 1 to 3) 
are equal to 1.0, 1.0 and 0.0 respectively. These assumptions 
are based on the fact that the clusters with a mass-loss of 
about 99% are mainly composed of SG stars and have ex¬ 
pansion factors less than 1.0, as explained in Section 3. We 
would like to stress that we interpolate in a 3D parameter 
space (? 7 , T, n/ai), so although all clusters are in the same 
tidal field they don’t have the same rtjai. 

In our analysis, we also consider the late-time adiabatic 
expansion of clusters as a result of the remnant gas expul¬ 
sion (for clusters with r > 10 ^) and also stellar evolution 
induced mass-loss by scaling the radii of all clusters accord¬ 
ing to the mass-radius relation of Hills (1980) which states 
that the radius of a cluster inversely scales with its mass, 
i.e. Ph oc assuming that the cluster remains in virial 

equilibrium after the initial significant mass-loss. We have 
used the AMUSE^ code (Portegies Zwart et al. 2009; Pelu- 
pessy et al. 2013; Portegies Zwart et al. 2013) and analytic 
stellar evolution models from Hurley, Pols & Tout (2000) 
to find the stellar evolution induced mass-loss for each clus¬ 
ter, which is on average equal to ~ 30% of the initial mass 
of the cluster, calculated for interval tend < t < 13.8 Gyr, 
where tend is the age of the cluster in physical units at the 
end of iV-body simulation. The IMF of EG and SG stars is 
a Kroupa (2001) IMF which extends to 100 M© and 8 M 0 
respectively. SG stars cannot be more massive than 8 M 0 , 
otherwise they will explode as SN and in the AGB scenario 
this will change the iron abundance of SG stars which is not 
consistent with observations (D’Ercole et al. 2008). 


^ AMUSE (Astrophysical Multipurpose Software Environment) 
is available at http://amusecode.org 


We use a least-squares method, to find the model which 
best matches the observations. We consider the distribution 
of Galactic GGs from the latest version of Harris (1996) cat¬ 
alogue in a 2D plane of mass vs radius and split this 2D 
plane into bins and calculate the normalized frequency of 
GGs in each bin, i.e. the number of GGs that are in each 
bin divided by the total number of GGs, so that we have 
a 2 D matrix of these normalized frequencies (hereafter O). 
We make a similar matrix for our simulated clusters (here¬ 
after S ). We then calculate the sum of the squared residuals 
between matrices O and S 

D = - Siif (10) 

ij 

By minimizing D, we can find the set of initial param¬ 
eters, given in Table 2, which best match the observational 
distribution. As an additional criterion, we only consider 
those set of initial parameters for which the distribution of 
the fraction of SG stars has a sample mean value of 50 ±5%. 
This way we can make sure that the fraction of SG stars in 
our simulated clusters are consistent with what we see in 
the observed clusters (D’Antona & Galoi 2008). In order to 
reduce the statistical errors, we do our MG simulations in 
the neighbourhood of each best-fitting model in the param¬ 
eter space for 20 different random seed numbers. We then 
take the mean values of D and the fraction of SG stars as 
the selection criteria. 

Table 3 lists our best-fitting models for different tidal 
fields, concentration of SG stars and dependence of clus¬ 
ter initial radii on the cluster masses. Fig. 4 illustrates the 
outcome of our MG simulations for the following models: 
(MD-RG2.0, A = 0.1), (MD-RG4.0, A = 0 . 2 ), (MD-RG8.5, 
A = 0.1) and (MI-RG8.5, A = 0 . 2 ) and compares them with 
the distribution of the observed clusters. As one can see, our 
best-fitting models match the distribution of the observed 
clusters very well. The mass and half-mass radius distribu¬ 
tions of our best-fitting models have roughly preserved their 
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logi„(M/Mo) logi„(rJpc) NJ(N, +N^) 


Figure 4. Comparison of the distribution of cluster masses (left), half-mass radii (middle) and the fraction of SG stars (right) for 
observed clusters (hatch-filled) and four of our best-fitting models listed in Table 3. Cluster masses and radii are taken from the most 
recent version of Harris (1996) and fraction of SG stars are taken from D’Antona & Caloi (2008). 



logio(M/M©) logio(r/, /pc) ^og^oiry /pc) 


Figure 5. Distribution of mass and half-mass radii of MD-RG8.5 (left and middle) and MI-RG8.5 (right) both with A = 0.2. We 
have plotted the best-fitting models when log(Texp/yr) is fixed and equal to 4.80 (dashed-blue line), 5.00(dashed-black line) and 5.20 
(solid-cyan line). The solid-green line shows the original best-fitting model when the gas expulsion time-scale is also a free parameter. 
The original values for log(Texp/yr) is 4.15 and 3.40 for MD-RG8.5 and MI-RG8.5 respectively. 


initial log-normal distributions and the mean of the distri¬ 
butions have shifted towards lower masses and higher radii 
respectively which is a direct consequence of gas expulsion. 
Due to the low number of observed GCs with measured MSP 
ratios, in our analysis we only fit the mean of the distribution 
of SG fraction and not the actual shape of the distribution. 
Fig. 4 shows that our best-fitting models have sample means 
of 50 zb 5% for the fraction of SG stars, which is consistent 
with observations. 

According to Table 3, the initial stellar masses of the 
GCs need to be of order 5 — 15 x 10^ M© with a gas frac¬ 
tion of at least r] = 1.0, meaning that for the significant 
mass-loss scenario to work we need as much mass in gas as 
in FG stars. In addition, GCs in stronger tidal fields need 
to be initially more massive compared to GCs in weaker 
tidal fields since mass-loss is stronger in stronger tidal fields. 
Eq. (8) gives an initial half-mass radius of ~ 1.0 pc for our 
best mass-dependent models roughly equal to that of mass- 
independent models. The required gas expulsion time-scales 
are all extremely short, Texp ~ 10^ yr (Tcross ~ 10^yr,T = 
0.1). To see if higher gas expulsion time-scales also lead to 
acceptable fits, we fixed the value of log (Texp / yr) to three 
different values of 4.80, 5.00 and 5.20 respectively and de¬ 
termine the best fitting-models for the MD-RG8.5 and MI- 
RG8.5 models with A = 0.20. Fig. 5 shows that for gas 
expulsion time-scales larger than T > 10^ yr, the final prop¬ 
erties of clusters are in strong disagreement with the ob¬ 
served clusters, especially for mass-dependent models, im- 



Figure 6. Fraction of SG stars as a function of the cluster mass 
for the simulated and the observed clusters (filled red squares). 
Simulated data points show an anti-correlation which is not seen 
in the observed clusters. In this plot, we show only a fraction of 
the simulated data points for clarity. 


plying that the gas expulsion time-scale must have been less 
than T = 10® yr. 

Our MG simulations also predict an anti-correlation be¬ 
tween the fraction of SG stars and the final mass of GCs as 
illustrated in Fig. 6. This is due to fact that to increase the 
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fraction of SG stars, GCs need to lose many of their FG stars 
and since FG stars constitute ~ 90% of the cluster initial 
mass, such GGs will have a lower final mass on average. We 
have used the Pearson correlation coefficient to quantify this 
anti-correlation 

_ ((a:-(a:))(2/-(2/))) 


where x and y corresponds to logio(^*/ ^o) ^ 2 /(W + 

N 2 ) respectively. The penultimate column of Table 3 shows 
the value of the anti-correlation for different models. As one 
can see this anti-correlation exists in all models, regardless 
of the tidal field strength or concentration of SG stars. The 
anti-correlation is more pronounced for mass-independent 
models with an average r value of ~ —0.85 compared to 
mass-dependent models with r ~ —0.68. Observed clusters, 
denoted by the filled red squares in Fig 6, only exhibit a 
weak anti-correlation with r = —0.07 as the fraction of SG 
stars is almost independent of the cluster mass. The 90% 
confidence interval of r for observed GGs, obtained from the 
bias-corrected and accelerated (BGa) bootstrap method of 
Efron (1987), is [-0.68, 0.61]. As a result the lower confidence 
bound for the anti-correlation coefficient of observed GGs is 
marginally in agreement with the mass-dependent models 
but it is still different from mass-independent models. 

This discrepancy can be explained in two ways. Either 
the significant mass-loss scenario does not work or we need 
to conduct more surveys to find the fraction of different stel¬ 
lar populations for more clusters. In either case, the exis¬ 
tence of such an anti-correlation could be used as a diagnos¬ 
tic to test our scenario. 

The analysis of Decressin et al. (2010) on the Baum¬ 
gardt k, Kroupa (2007) models show a similar relation be¬ 
tween the fraction of SG stars and the number of bound 
stars. Since the number of bound stars is proportional to the 
total mass of the cluster, the result of their work matches 
the anti-correlation that we see in Eig. 6. 


6 CONCLUSIONS AND DISCUSSION 

Using a large grid of A-body and Monte Carlo simulations 
we have studied the consequences of primordial mass-loss for 
GCs with MSPs to put constraints on their initial conditions 
and find the best match with observations. We have demon¬ 
strated that primordial mass-loss is able to simultaneously 
reproduce the present-day distribution of GCs in the mass- 
radius plane (Eig. 4) and explain the large fraction of second 
generation stars in the cluster. However this is only possi¬ 
ble if: (1) the total mass of the gas remaining in the cluster 
is equal to that of first generation stars M ~ 10® M© and 
(2) very short gas expulsion time-scales of less than 10® yr 
which is equal to about one initial crossing time. In this 
case, typical initial masses of globular clusters are around 
M = 10® M© and their initial half-mass radii are around 
r/i = 1 pc (Table 3). 

According to Decressin et al. (2010) for a gas cloud with 
an initial mass of ~ 10® M© and an initial half-mass radius of 
0.5 pc, around 50 SNe are needed to unbind the gas cloud. 
In our case the gas clouds, which are more concentrated 
than the first generation stars, have initial half-mass radii 
of 0.1 and 0.2 pc. Since the gravitational potential energy 


scales with oc R~^ ^ our gas clouds will need at least of order 
125 — 250 SN explosions to be dispersed. According to Eig. 

5 of Decressin et al. (2010) in a cluster with M = 10® M©, 
at most 400 SNe/Myr will explode within 1 Myr implying 
that only 40 SNe will go off in 10® yr which is a factor 3 to 

6 below the required limit. As a result SN explosions seem 
not to able to generate enough energy to expel the gas over 
the short time-scales we need in our scenario. 

In addition, Krause et al. (2012) have shown that the 
Rayleigh-Taylor (Sharp 1984) instability destroys the huge 
gas shells (superbubbles) made by SN explosions before they 
build up enough speed to leave the cluster, thus such super¬ 
bubbles are ineffective in expelling the gas even if they have 
enough energy. Instead Krause et al. (2012) propose accre¬ 
tion onto dark remnants, such as neutron stars and black 
holes, as a promising mechanism which is capable to over¬ 
come the Rayleigh-Taylor instability and lead to rapid gas 
expulsion. According to their model, the dark remnants be¬ 
come active after the SN phase (t > 35 Myr) (Krause et 
al. 2013) and are able to unbind the intra-cluster medium 
in 0.03 Myr to 0.06 Myr depending on if neutron stars also 
contribute to the gas expulsion (Krause et al. 2012, 2013). 

According to Krause et al. (2012, 2013), This model 
works for protoclusters whose masses are less than 2 x 
10^ M© above which the gas cannot be ejected and will be re¬ 
tained in the cluster. This is intriguing, because firstly such a 
mass limit encompasses the majority of GGs except for very 
massive ones such as u Gen, which is 10 times more massive, 
and secondly MSPs in such massive GGs have different iron 
abundances which could be a result of gas retention (Krause 
et al. 2012). 

The gas expulsion time-scales and the initial mass of 
the gas clouds the we obtain in this work, match results by 
Krause et al. (2012, 2013) very well. However it is not clear 
whether GGs can retain gas clouds with masses of order 
~ 10® M© for ~ 35 Myr. Observations of young massive star 
clusters by Bastian, Hollyhead k Gabrera-Ziri (2014) and 
Hollyhead et al. (2015) show that they have cleared out their 
natal gas within a few Myrs. If it was the case for GGs 
as well, then it poses a serious challenge for the scenario 
proposed by Krause et al. (2012, 2013) and would imply 
that either stellar winds and supernova explosions do have 
to expel the gas or that the gas is completely consumed into 
stars after a few Myr in which case the scenario suggested 
here would not work. 

Another possibility is that only the centres of star clus¬ 
ters are gas free, but that the gas is still present in the outer 
regions. If converted to a physical scale, the half-mass radius 
of the second generation stars in our best-fitting models is 
around ~ 0.1 pc which is still far less than half-mass radius 
of the whole cluster (~ Ipc). As a result, it is possible to 
start with clusters that have a high star formation efficiency 
and little gas in the very centre, followed by a region with 
considerable amounts of remaining gas at intermediate radii 
and then the first generation stars at large radii, and still 
end up with large numbers of second generation stars after 
gas expulsion. This should work as long as first and second 
generation stars are well separated in space, i.e. a 2 /ai <C 1. 

In the AGB scenario, second generation stars form after 
SN explosions or dark remnants have expelled the primordial 
gas not accreted into stars from the clusters. As a result, the 
clusters need to accrete significant amounts of unprocessed 
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new gas into their centers to start formation of second gen¬ 
eration stars, while at the same time preventing the dark 
remnants to immediately eject this gas. According to our 
scenario, when it finally happens, the gas expulsion has to 
be very rapid. At the moment it is completely unclear if 
and how this is possible. In addition, accretion of gas is only 
possible for clusters which move with a low relative veloc¬ 
ity to an surrounding gaseous medium (Pflamm-Altenburg 
& Kroupa 2006), however second generation stars have been 
found in almost all massive globular clusters, independent of 
their orbits and position in the Milky Way. These problems 
do not exist in the FRMS scenario or any other scenario 
that form second generation stars within a few Myrs of the 
first generation ones, since the gas out of which the second 
generation forms is already in the cluster. 

Since the fraction of second generation stars increases 
as a result of long-term dynamical evolution of the clusters 
in the galactic tidal field (Decressin, Baumgardt & Kroupa 
2008), the fraction of second generation stars at the end of 
the gas expulsion phase could be lower than the present-day 
fraction. This would mean that the gas expulsion time-scales 
could be larger than what we found here since the fraction 
of second generation stars at the end of the gas expulsion 
decreases with the gas expulsion time-scale. However since 
the lifetimes of most globular clusters are significantly longer 
than a Hubble time (Baumgardt & Makino 2003), we do 
not expect the fraction of second generation stars to change 
significantly over a Hubble time due to dynamical evolution, 
and therefore our upper limit of 10^ yr for the gas expulsion 
time-scale is unlikely to change significantly. 

The outcome of our simulations shows that fraction of 
second generation stars is inversely proportional to the fi¬ 
nal cluster mass (Fig. 6). This anti-correlation, which is in 
agreement with Decressin et al. (2010), is one of the impli¬ 
cations of the primordial mass-loss and can be used to test 
the feasibility of this scenario. Observations show that such 
an anti-correlation, albeit weaker, also exist in the Galactic 
GCs. For our simulated clusters the anti-correlation coef¬ 
ficient ranges from —0.50 to —0.87, whereas for Galactic 
GCs it is about r ~ —0.07 with a 90% confidence interval 
of [—0.68,0.61]. As a result Galactic GCs show a relatively 
weaker anti-correlation. However, given the 90% confidence 
interval on the correlation coefficient of Galactic GCs, the 
data is also consistent with no anti-correlation or even a 
positive correlation. 

The discrepancy between theory and observation might 
be due to low-number statistics. Also in individual GCs, the 
number ratio of second generation stars has been measured 
only over a limited range in radius and observations show 
that the ratio is varying with radius (Lardo et al. 2011). 
Better observations are therefore needed to test if an anti¬ 
correlation similar to the one predicted by our models exists 
in Galactic GCs. 

In our Monte Carlo simulations we have assumed a log¬ 
normal distribution for the initial cluster mass function. In¬ 
stead of a log-normal relation, one can also assume a power- 
law distribution dA^ oc M~^ for the cluster mass function, 
with a ^ 2, as seen for young massive star clusters in inter¬ 
acting and merging galaxies (Whitmore & Schweizer 1995; 
Whitmore et al. 1999). Baumgardt, Kroupa & Parmentier 
(2008) for example studied the effect of residual gas expul¬ 
sion on gas embedded star clusters and found that it is pos¬ 


sible to turn a log-normal mass function into a power-law 
over a Hubble time due to gas expulsion. They found that 
this effect is almost independent of the strength of the ex¬ 
ternal tidal field or the assumed model for gas expulsion. As 
a result, our model could also work for an initial power-law 
distribution. A potential problem for such a mass function 
could be the over-production of the field halo stars due to 
the large number of disrupted clusters. A detailed numeri¬ 
cal analysis of this effect is beyond the scope of the present 
paper and can be the subject of a future paper. 
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